;; Copyright 2009 Adam C. Foltzer

;; Licensed under the Apache License, Version 2.0 (the "License");
;; you may not use this file except in compliance with the License.
;; You may obtain a copy of the License at

;; http://www.apache.org/licenses/LICENSE-2.0

;; Unless required by applicable law or agreed to in writing, software
;; distributed under the License is distributed on an "AS IS" BASIS,
;; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
;; See the License for the specific language governing permissions and
;; limitations under the License.

(define mm-entrywise-op
  (lambda (op)
    (lambda m*
      (let ([map-fn (lambda i*
                      (if (vector? (car i*))
                          (apply mm-matrix-sum i*)
                          (apply op i*)))])
        (apply vector-map map-fn m*)))))

(define mm-matrix-sum
  (mm-entrywise-op +))

(define mm-matrix-diff
  (mm-entrywise-op -))

(define mm-entrywise-product
  (mm-entrywise-op *))

(define mm-dot-product
  (lambda (v1 v2)
    (let ([v-length (vector-length v1)])
      (let loop ([index 0] [acc 0])
        (if (= index v-length)
            acc
            (let* ([v1-elt (vector-ref v1 index)]
                   [v2-elt (vector-ref v2 index)]
                   [prod (* v1-elt v2-elt)])
              (loop (add1 index) (+ acc prod))))))))

(define mm-magnitude
  (lambda (v)
    (sqrt (mm-dot-product v v))))

(define mm-scalar-multiply
  (lambda (m a)
    (let ([map-fn (lambda (i)
                    (if (vector? i)
                        (mm-scalar-multiply i a)
                        (* i a)))])
      (vector-map map-fn m))))

(define mm-normalize
  (lambda (v)
    (mm-scalar-multiply v (/ (mm-magnitude v)))))

;;Takes two functions f1 f2 and an arbitrary number of vectors of equal length.
;;From index zero to the end of the vectors, applies f1 to the vector elements
;;at that index, then applies f2 to the result and the result of the recursion.
(define vector-zip-fold
  (lambda (f1 f2 . v*)
    (let ([v-stop (sub1 (vector-length (car v*)))])
      (let loop ([index 0])
        (let* ([this-row (map (lambda (v) (vector-ref v index)) v*)]
               [this-result (apply f1 this-row)])
          (if (= index v-stop)
              this-result
              (f2 this-result (loop (add1 index)))))))))

;;Euclidean distance between two points (represented by vectors)
(define mm-distance
  (lambda (p1 p2)
    (let* ([map-fn (lambda (p q) (expt (- p q) 2))]
           [sum (vector-zip-fold map-fn + p1 p2)])
      (sqrt sum))))

